From individual to collective behaviour of coupled velocity 
jump processes: a locust example 

Radek ErbaiJ^ Jan HaskovecH 

Abstract. A class of stochastic individual-based models, written in terms of coupled velocity jump 
processes, is presented and analysed. This modelling approach incorporates recent experimental findings 
, on behaviour of locusts. It exhibits nontrivial dynamics with a "phase change" behaviour and recovers 
' the observed group directional switching. Estimates of the expected switching times, in terms of number 
O • of individuals and values of the model coefficients, are obtained using the corresponding Fokker-Planck 
. equation. In the limit of large populations, a system of two kinetic equations with nonlocal and nonlinear 
[ right hand side is derived and analyzed. The existence of its solutions is proven and the system's long-time 
behaviour is investigated. Finally, a first step towards the mean field limit of topological interactions is 
made by studying the effect of shrinking the interaction radius in the individual-based model when the 
I number of individuals grows. 
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1 Introduction 

K*" ■ Individual-based behaviour in biology can be often modelled as a velocity jump process [20]. Here, the 
^ , velocity of an individual is subject to sudden changes ( "jumps" ) at random instants. If the velocity changes 
I are completely random, this process simply leads to diffusive spreading of individuals in an appropriate 
' limit [18]. The situation becomes more complicated whenever the velocity changes are biased according to 
'si" . an individual's environment. A classical example is bacterial chemotaxis 0. Individual bacteria change 
their frequency of velocity changes according to their environment. If they swim in a favourable direction 
(e.g. towards a nutrient source), they are less likely to change their direction. On the other hand, they 
are more likely to turn if they are heading away from a foodstuff [9] . 

In this paper, we modify the velocity jump methodology to model the behaviour of locusts. Our model 
is motivated by the recent experiments of Buhl et al [2] . They studied an experimental setting, in which 
■ locust nymphs marched in a ring-shaped arena. The collective behaviour depended strongly on locust 
density. At low densities, there was a low incidence of alignment among individuals. Intermediate densities 
were characterized by long periods of collective motion in one direction along the arena interrupted by 
rapid changes of group direction. If the density of locusts was further increased, the group quickly adopted 
a common and persistent rotational direction. Yates et al [21] analysed experimental data of Buhl et al [2] 
and proposed that the frequency of random changes in the direction of an individual increases when the 
individual looses the alignment with the rest of the group. In this paper, we incorporate this observation 
into a stochastic individual-based model formulated as a velocity jump process. We show that this model, 
although phenomenologically very simple, has the same predictive power as other modelling approaches 
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previously used in this area [3 [2]. In particular, it exhibits (i) a rapid transition from disordered movement 
of individuals to highly aligned collective motion as the size of the group grows, and (ii) sudden and rapid 
switching of the group direction, with frequency decreasing as group size increases. 

The individual-based model is introduced in Section [2l The ring-shaped arena, used in Buhl's experi- 
ments [2], is modelled as one-dimensional interval with periodic boundary conditions. Locusts march with 
a constant speed and each individual switches its direction randomly. The individual switching frequency 
increases in response to a loss of alignment. In Section [Sj the corresponding Fokker-Planck equation is 
derived for the system with global interactions and possible types of qualitative behaviour of the system 
are classified. For the case of ordered group motion, where two distinct metastable states exist, an ap- 
proximate analytic formula for the mean switching time between these two states is derived. Then, in 
Section HJ the kinetic formulation of the model is obtained in the limit as the number of locusts tends to 
infinity. The existence of solutions of the kinetic model is shown in Section [5] and the long time behaviour 
is investigated in Section [6l We conclude with analysis of the dependence of collective behaviour on the 
size of the interaction radius of the individuals in Section [71 



2 Individual based model 

We consider a group of N agents (locusts) with time-dependent positions Xi{t) and velocities Vi{t), i = 
1, . . . ,N. To mimic the ring-shaped arena set-up of [2], we assume that the agents move along a one- 
dimensional circle, which we identify with the interval = [0, 1) with periodic boundary conditions, and 
move either to the right or to the left with the same unit speed, i.e. 

Xiit) Vi£ {-1, 1} and -^{t) = Viit). (2.1) 

We define the local average velocity of the ensemble, seen by the i-th. agent, as 

^loc ^ J2m=l "^il^i - Xm\)Vm 2) 
Yll=lW{\xi-Xm\) 

where tt; is a weight function defined on 0, with the properties: 
[Al] w is bounded and nonnegative on 0, 
[A2] w{0) > 0. 

For example, w(s) = xp.o-jCs), where X[o,a] is the characteristic function of the interval [0,a] and a > 
is a interaction radius, satisfies conditions |[A1]| and |[A2][ This is a common choice of w in biological 
applications [71[23]. It is worth noting that, due to the assumption | [ A2] | the definition ()2.2p always makes 
sense and ul""^ G [—1)1]- 

The agents switch their velocities to the opposite direction (i.e., from Vi = 1 to Vi = —1 and vice 
versa) based on N independent Poisson processes with the rates 

li = lo + bavi-u'f''), i = l,...,N, 

where 70 > and 6 > are fixed parameters and the "response to disalignment" function ^ : [—2, 2] — ?> 
[0, 00) is assumed to be convex, differentiable and symmetric with respect to the origin. Taking the Taylor 
expansion of ^(s) around s = 0, we obtain 

C(s) = ao + a2S^ + 0{s^). 
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N=5 N=7 N=12 




Figure 1: An example of the transition to ordered motion as N grows. We use l\2.1\i - (K4^ with 6 = 1, 
7o = 0.2 and w = X[o,o.2] ■ Shown are the normalized histograms of the group mean velocities u = jj "^iLi Vi 
recorded in 10^ time steps of length 10~^ , with N = 5 (left panel), N = 7 (middle panel) and N = 12 
individuals (right panel). The system does not prefer any particular state for N = 5. Two quasi-stable 
states of ordered collective motion are easily recognizable for N = 12. 

We can set oq = without loss of generality, because it can be absorbed in 70. Since the individuals 
switch their velocities less frequently when they are aligned ([21]), ^(-s) has a global minimum at s = 0. 
This implies that 02 > 0. If 02 > 0, we can set a2 = 1 by choosing an appropriate time scale. Therefore, 
^(s) has the general form + 0{s^). For the rest of the paper, we choose the form ^(s) = for 
simplicity. Other choices are certainly possible, for example, in the limiting case 02 = the leading order 
approximation is given by a higher order term, which, however, complicates the analysis. However, it is 
worth noting that the derivation of the kinetic equation performed in Section H] is possible, for example, 
also for ^(s) = \s\'^, n > 3. 

With ^(s) = s^, the turning rate "from the right to the left", J^^^, and the rate for the opposite 
turn, ryL^R^ a,re given by 

7f = 70 + 6 (1 - u\°''f for the switch from = 1 to = -1 , (2.3) 
^f'^^ = ^0 + 6(1 + u\°''f for the switch from v-i = -1 to Vi = I . (2.4) 

This velocity jump process describes the tendency of the individuals to align their velocities to the 
average velocity of their neighbors. Despite its relative simplicity, the model provides similar predictions 
as the Vicsek and Czirok model [7] and its modification [23] and is in qualitative agreement with the 
experimental observations made in [2] : the transition to ordered motion as N grows (Figure [T]) and the 
density-dependent switching behaviour between the ordered states (Figure [2l bottom). 
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Figure 2: An example of the behaviour of the model l\2.1\i - ^274\ i with global interactions (w = 1): the large 
noise case (top) with N = 20, b = 1 and jq = 1.3 and small noise case (bottom) with N = 20, 6=1 and 
7o = 0.3. Left are the histograms of the group mean velocities u recorded in 10^ time steps of length 10~^, 
compared to the plot of the (properly scaled) stationary solution ps of the corresponding Fokker-Planck 
equation (solid line). Right are the plots of the temporal evolution of the group mean velocity u during 
5 X 10^ timesteps. In the small noise case, one can clearly distinguish the two quasi- stationary states and 
observe the switching between them. 
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3 Analysis of the individual based model with global interactions 

In this section, we simplify the individual-based model by assuming w = 1 in (1221), i.e. u'f" = u for all 
i = 1, . . . , N, where 

1 A , , 2r(t) -N , , 

4 = 1 

where r(t) is the number of individuals which are going to the right at time t. Using this simplification, 
we will obtain an explicit formula for the mean switching time between the ordered states. However, 
for the derivation of the kinetic decription and its analysis (Section U]), we will allow general weights w, 
imposing only the assumption [ Al ] and a slightly reinforced version of [A2][ 

Let p{r,t) be the probability that r individuals move to the right (i.e., with velocity 1) at time t>0. 
It satisfies the master equation 

^p(r,t) = (r + l)7«^^p(r + l,t)-r7^^V^t) 

+ {N -r+l)-y^~^^p{r -l,t) - (N -r)j^-^^p{r,t) , (3.2) 

where j^^^ and 7^^^ are given by ()2.3p and (12. 4p . respectively. The subscript i in ()2.3p - (l2.4p is dropped 
in (|3.2p because all individuals have the same turning rates. Using the system size expansion [22 and the 
definition (|3.ip of the average velocity u{t), we obtain the following Fokker-Planck equation 



dpiu, t)_d_ _ _ ^^^^ ^^l^2_ ^ _ ^^^^ 



dt 



(3.3) 



where p{u, t) is the probability distribution function of the average velocity p.ip at time t. The stationary 
solution ps of (|3.3p is 

Ps{u) = Cexp[-<^Niu)] (3.4) 
where C is the normalization constant and the potential <I>Ar is given by 

^n{u) = -^u2 + (1 - ^) In (70 + 6(1 - n2)) . (3.5) 

The comparision of the stationary probability distribution function ps with the results obtained by long- 
time simulation of the stochastic individual-based model is shown in Figure [21 Differentiating (j3.5p twice, 
we obtain 

Consequently, we distinguish the following two cases: 
7q 2 

(1) Large noise: If > 1 + — , then has the global minimum at u = 0. 

(2) Small noise: If ^ ^ "f" ' ^^^^ ^ local maximum at u = 0. The only local and global 
minima are at it n, where 



2 70 
IH -. 

N b 
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It is worth noting that > 1 if is smah, namely for < 26/79. This is a consequence of approximations 
made during the derivation of the Fokker-Planck equation (13. 3p . 

In the large noise case, the system prefers the disordered state u = 0, while, in the small noise case, 
the system has two preferred ordered states ±n<j where it spends most of the time and eventually switches 
between them (see Figure [21 bottom panel). 

Using Kramers theory \15\ I13j . the mean switching time tn between the states —Ug and Us can be 
approximated as 

( ^ , Ntt exp($jv(0) - ^jv(-n,)) 
Ta[(—Us i-> Us) ~ r = — , 3.7) 

which has an exponential asymptotic with respect to large N given by 



b--foy 2(70 + b) 



&-70 _ To 270 



70 + 6 



This is in agreement with the experimental observations, [2], as well as with the modified Czirok-Vicsek 
model of [23], where the mean switching time is as well exponential in A^. Finally, it is interesting to note 
that with the transform b = 796 in ()3.7p . one has 

r(iV,7o,6) =7or(A^,6), 

i.e., if 6/70 is kept fixed, the mean switching time scales linearly with 70. 

In the numerical experiment shown in Figure [2] bottom (small noise case with 70 = 0.3, 6 = 1 and 
N = 20 agents), we have two metastable states located approximately at Us = ±0.894. The estimate 
mean turning time given by formula (|3.7p is tn{—Us ^ Us) = 61.1. Performing 10^ time steps of length 
10~^, the observed mean switching time (defined as a mean transition time between the states v = —0.8 
and V = 0.8 or vice versa) was 58.8, showing a very good agreement. 



4 Kinetic description 

In this section we derive the kinetic description of the system of A'^ interacting agents and formally pass 
to the limit A^ — ?> 00 to obtain the corresponding kinetic equation. For this, we have to accept a slight 
reinforcement of the assumption A2]on w, namely, 

[A2'] w > on the interval [0, r) for some r > 0. 

The state of the system of A^ agents at time t > is described by the probability density function 
{t,xi,vi, . . . ,X]\f,VN) of finding the i-th agent in position G with velocity Vi G {—1,1}, for i = 
1, . . . ,N. When convenient, we will use the abbreviation = p^{t, x, v), with x = (xi, . . . ,xn) G 
and V = [vi, . . . ,vn) G {±1}^. The probability density p^(t,x, v) evolves according to 

|^'^ + E-^^^'^ = -^[^'^] + ^b'^]' (4-1) 

1=1 * 

where /^[p"^] (resp. ^[p^]) is the loss (resp. gain) term corresponding to the velocity jumps. Using 
()2.3p - ()2.4p . the rate of the switch Vi —Vi is for z G {1, ... , A^} given by 

7i(x, v) = 70 + 6 (tij - u[xi;x,v])'^ , 



6 



where u[2;j;x, v] = v}°^ is defined by (j2.2p . i.e. u[2;j;x, v] is the local average velocity seen by an agent 
located at Xi € $7, based on the system configuration [x, v]. Consequently, the loss term is given by 



TV 



/:[p^](f,x,v) = J]7,(x,v)p^(t,x,v). (4.2) 

i=l 

We define the operator Mi : {±1}*= {±1}'' for A; € {1, ... , N} and i G {1, . . . , k} by 

Mi(v) = {vi, . . . ,Vi-i, -Vi,Vi-i, . . .,Vk), (4.3) 

i.e. Mj(v) denotes the velocity vector created from v by changing the sign of its i-th. component. Then 
the gain term Q[p^]{t,'x.,v) is given by 

TV 

g[p^](t,x,v) = J]7.(x,M,(v))j>^(t,x,M,(v)), (4.4) 

i=l 

where 

7j(x,Mj(v)) = 70 + & (- - u[xi;x, Mi(v)])^ . 
Consequently, the right hand side of (|4.ip is 

(g-/:)[p^](x,v) = 7o(^J^p^(x,Af,(v))-iVp^(x,v)^ (4.5) 

TV 
1=1 

where we dropped the dependance on time t to simplify the notation. Finally, we postulate the so-called 
indistinguishability-of-particles: We only consider solutions p^ that are indifferent to permutations of their 
(xj, Uj)-arguments. Such solutions are admissible, since the equation (j4.ip with the collision operator ()4.5p 
is as well indifferent with respect to interchange of the (xj, 'L'j)-pairs. 

4.1 Derivation of the BBGKY hierarchy 

To derive an analogue of what is called the BBGKY hierarchy in the classical kinetic theory of gases (see, 
for instance, [6]), we define, for A; = 1, . . . , A^, the fc-agent marginals 

pN,k^^k^^k^^ y f ^N^^k^-^^k-^^-^ x'^ G 0^ G {±l}^ (4.6) 

VG{±1}^~'= 

In what follows, the coordinates of x'^ and x will be denoted as 

x'^ = {Xi, X2, . . . , X^) , X = (xi, X2, . . . , Xj^-k) , 
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that is, X = (x^,x) = (x^;', x^, • • • , a;^, xi, X2, . . . ,^Ar_fc). The same notational convention is used for 
velocities, i.e., v = (v'^', v) = {vi,V2, ■ ■ ■ ,v^,vi,V2, ■ ■ ■ ,VN~k)- Note that, due to the indistinguishabihty- 
of-particles, the marginals are well defined and indifferent with respect to permutations of the pairs of 
arguments {xi,Vi), i = 1, . . . , k. Integrating (j4.ip . we obtain 

|p^''=(x^v'=)+^.,Ap^'.'=(x^v'=)= L /a-^)[p^](x^x,v^v)dx. (4.7) 

Substituting (j4.5p into the right-hand side of (j4.7p . we get 

I -, AT I. 

k 
i=l 

+ ^ E ^E{(-^«'-^[^-^^^'^^*(^')'^])V(x^x,M,(v^),v) 

— ^ r \ t ^ N —I- ^ ^ A — 1 



ve{±i} 



- (uf - ti[xj;x'^,x, v'^, v])^p^(x'^,x, v*^, v)| 

k 

= (70 + 6) ^(p^'^-(x',M,(v'=)) -p^'^(x^v^)) 
1=1 

+ 26 E / V«K^x[xi;x^x,M,(v'=),v]p^(x^x,Mi(v^),v) 

+ ti[xi;x'^,x, v'^^ v]p^(x''\x, v'^, v)^ dx 
+ ^ E ^E(H^-^''^'^^^(^')'^])'^''^(^''X,M^(v'^),v) 



- (n[x,;x^x,v^v])'p^(x^x,v^v)j dx. (4.8) 

Since we are interested in the limit N oo with k fixed, we can rewrite the definition ()2.2p as follows 

k ^ k ^-^ _ J2t=l'^(\Xm-Xi\)vt + ElZlW{\Xm-Xi\)Vm 



u\Xi;x ,x, V ,vj — ^ j^—j^ — 

Em=l ^(l^^m + Ylm=l ^(l^m " Xi\) 

-- n[xi;x,v] + 0( 41 , (4.9) 



N 



where we define 

yN-k 



n[z;x,v]= ^-^"^'"--"l^^- for.EO. 

22ra=lWi\^rn- Z\) 
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Substituting (gS]) in (jilSj) and (jO]) . we obtain the BBGKY hierarchy 



^pN,k^^k ^Y^^.dN,k^^k = (70 + 6) ^(p^''=(x^M,(v'=)) -p^''=(x^v* 

L/C _ CJJi"i ^ ^ 



i=l ' i=l 



+ 2h^v^^ (q^'^(xi] x^Mi(v'=)j x^v'=jj (4.10) 

1=1 

+ (^^; x^M,(v'=)) - r^'^(x.; x^ v'^)) + O (A) , 



where 



ve{±i} 



I 'u[z;x,v]p^(x^x,v^v)dx. (4.11) 



and 



„N,k 



z; x'^jV^j = / 'u[z;x, v]^p^(x^',x, v'', v) dx. (4.12) 



^e{±i} 



4.2 Passage to the limit N ^ oo 

The usual procedure of deriving the mean field equation is to write the BBGKY hierarchy (j4.10p in terms 

Qf pN,k g^j-j^ 

pass to the limit N — y oo to obtain the so-called Boltzrnann hierarchy for : — lim^—^.ooP^'''^ 
|22j . Then, one shows that the Boltzmann hierarchy admits solutions generated by the molecular chaos 
ansatz (see below). In our case, however, this strategy cannot be pursued; although we could derive 
uniform estimates allowing us to pass to the limit — )• cx) in the BBGKY hierarchy, we are not able to 
express the limiting 

^N,k terms of p^. Consequently, 
we have no clue what the correct molecular chaos ansatz for and r'^ should be. 

Instead, we assume the propagation of chaos already at the level of the BBGKY hierarchy, before 
passing to the limit N ^ oo: we assume that, for large N, is well approximated by the product of the 
limiting one-particle marginals p '. — liniTv^QQ p^'\ i.e. 

N 

p^{t,x,v) ^Ylp{t,Xi,Vi) for alH > 0, x e Jl^, V e {±1}^. (4.13) 

i=l 

This corresponds to vanishing statistical dependence (correlations) between the agents as ^ oo and 
is the usual phenomenon observed in systems of interacting particles, see for instance [B] in the context 
of classical kinetic theory or [T7] in the context of biological systems. Moreover, if one interprets 
qN,k (^j-ggp^ ^N,k^ ^l^g £j,g|. (^pggp^ second) order moment of p'^ with respect to tt[z;x, v], then one can 
understand (j4.13p as the moment closure assumption for the non-closed system of moments generated 
by gl}. 

The essential point is that now we may insert ()4.13p into (j4.1ip and (j4.12p to obtain explicit expressions 
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for and in terms of p\ 

g^(z;x^v'=)= lim q^^^iz'^^^ 

N-^oo \ 



lim q'^'^iz; x'',v'' = (4.14) 



and 



"^(z;x^v'=)= lim r^'^fz;x^v^' 

np(x^.f)') X lim 5: / f^fe%^%T np(-^,^odx. 



r'^(z;x^v'') = lim r^"''^ z; x^vM = (4.15) 



\±x'^ jv^oo ^ QN-k \ V^^"" infix - zl 

Study of the limit iV ^ oo in (lil^ and (lil^ 

We start by setting k = 1, which is the case considered in Section [4.31 To simplify the notation, we drop 
the bars over x and v, and, without loss of generality, choose z = 0. First, we explore the symmetry of 
the expression (j4.14p as follows: 

= yZ Tl^T — r^^—r^iPi^niA)-p{xm,-l)] V TT p(2;i,Vi)dx 



' Jf^iv-i (iV - l)57v(x) , , 

/ '^nc / J i^m) n £'(xi) dx 

Jn^-i [N - l)5Ar(x) .-^-L 

m=l / \ / j^m 



with the notation 



£»(x) := p(x, 1) +p(x, -1) , 
j(x) := p(x, 1) -p(x, -1) , 



and 



N-l 

Due to the normalization of p^, we have f^Q{x)dx = 1. Consequently, in what follows we denote by 
Pg(t) the time dependent probability measure corresponding to the probability density Q{t). Since w is 
bounded and nonnegative by assumption | [Al ] it is integrable with respect to and we may define 

/ := / w{x)dPi,{x) > 0. (4.17) 
Jn 
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The forthcoming analysis will be performed given the assumption / > 0; the case / = will be discussed 
in Remark [21 

Lemma 1 Let Pg he a probability measure on $7 with density g and j £ L^{Q) such that \j\ < p. Let 
w : Q —?■ [0, oo) with w € L°°{Q) be such that the integral I defined by ( |^. j7[ ) is positive. Define 

f w{xi)j{xi^ ^ ^ 



Qn '■- 



1=2 



Then 



lim Qn = j[ w{y)j{y)dy. 



Proof: We can treat w{y) as a random variable with respect to the probability measure Pg{y)- The 
essential tool of the proof is the law of large numbers, which states that Sn{^) converges to / in measure, 
in the sense that for each e > 0, 



lim -1 ({x € 0^-1; |5^(x) - I| > e}) = , 



(4.18) 



where P^ denotes the {N — l)-fold tensor product of the probability measures P^. Moreover, the 
existence of the m-th order moment of w with respect to Pg, 

\w{y)r dPgiy) < ^ , 

n 

implies the rate of convergence (see [l] ) 

lim {N - l)'"-ipj^-i ({x G n^-^; |5^(x) - I| > e}) = . 

Let us denote by A]^{e) the set {x € fi^"^, |5Ar(x) — /| < e} and by A'f^{e) its complement in fl 
Choosing < e < L/2 and x G 74jv(e), we have the estimate 



(4.19) 



w{xi) w{xi) 



Consequently, 



Ajv{e) 



5Ar(x 
w{xi) w{xi) 



2e , , 
< jiw{xi) . 



L 



dPf-l(x)<| / ^Xi)dPf"Hx) = y. 



On the other hand, the integral over A'^{e) is estimated with 



w{xi) w{xi) 



Sn{x.) 



I 



dPf-i(x) < 



w{xi) 



— dP-- 



(x)+ / 
J A 



w{xi) 

AUe) 'S'iv(x) 



dP 



N-l, 



X 



The first term on the right hand side converges to zero as A'^ — )• oo by (j4.18p and the boundedness of w. 
Using ()4.16p . the second term is estimated by 



JaUs) 5'7v(x) ^ 



\^)<{N-l)P^'\A%{e)) 
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and vanishes as ^ oo due to (j4.19p with m = 2. Consequently, we have shown that the term 

N-l 



Qn - J [ w{y)j{y)dy 
^ Jn 



< 



w{xi) w{xi) 
QN-i \Sn{x) I 

wi^xi) w{xi) 



j(a:;i)dxi JJ dPg 



1=2 



I 



dPf-Hx) 



can be made arbitrarily small, for sufficiently large A^. 



The formula (j4.15p can be analysed in a similar way as we did for (j4.14p . Namely, we exploit its symmetry 
as follows. 



E 



Jn^-r I (iV-l)5^(x) 



^2 N-~l 



Y[ pixi,Vi)dx 



VG{±1P 

(Af-1)2 ^ ^ J^N-l 



i=l 



w{Xi)w{Xm)jixi)jiXm) 



dxidxm JJ dPg{xk) 



+ 



- 1)2 



N-2 f w{xi)w{x2)jixi)jix2) 



N-l 



1 f w{xi)'^ 



dxi dx2 Y[ ^Pgi^k) 



k=3 



dPi^-i(x). 



N - 1 J^N-i S2,(x) 
The limiting behaviour of the first term is studied in the following Lemma: 
Lemma 2 With the assumptions and notation of Lemma [7], and defining 

w{xi)iv{x2)jixi)j{x2) 



N-l 



we have 



■d2;idx2 Yi dP^(a^fc) ; 



A:=3 



(4.20) 



lim i^jv = 7 / w{y)j{y)dy 
N^oo \I Jn 

Proof: We follow the lines of the proof of Lemma [TJ Defining again Aj\f{e) := {x € ft^"^, |5'Ar(x) — / | < 
e} and A^(e) := \ An{£), with < e < 1/2, we derive the estimates 



■w{xi)'w{x2) w{xi)'w{x2) 



p 



We 

T 
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and 



w{xi)w{x2) w{xi)w{x2) 



where the first term vanishes in the hmit N oo due to (j4.18p and the second one by (j4.19p with m = 3. 

■ 

By a slight modification of the above Lemma, one obtains the hmit of the second term of ()4.20p . namely 

lim / ^ dP^H^) = i / w{yf dPM . 



Therefore, the limit as iV oo of (j4.20p is 



lim I ^rTw'^rlT n P^^^^ = hm RN=(\j w{y)3{y) dy) , (4.21) 

where we used Lemma [2j 

Remark 1 Lemmas d and \M were formulated for the case k = 1, however, all the calculations can be 
easily generalized for any (fixed) value of k. 

4.3 Derivation of kinetic and hydrodynamic description 

Let us denote p~^{t,x) := p{t,x, 1) and p^{t,x) := p{t,x, —1). We define u{t,x), the continuous analogue 
of the local average velocity (j2.2p . by 

uit,x):=^^ iw , . ^ , for x^So\p^,p \{t), (4.22 

j^w{\x - z\){p+{t,z) +p {t,z))dz 

where 

Soy,p-]{t) := j^w{\x-z\){p+{t,z)+p-{t,z))dz = 0^ . 

We extend the definition of u{t, x) to the whole domain Q by setting u{t, x) = for x G So\p~^ , p~]{t) , see 
Remark [21 

The kinetic equation for p is obtained by setting = 1 in the BBGKY-hierarchy (j4.10p with the 
molecular chaos assumption (j4.13p and passing to the limit N ^ oo. Using Lemma [1] and formula (|4.2ip . 
we obtain 

q^{x;x,v) = p{t,x,v)u{t,x) , and r^{x;x,v) = p{t,x,v)u{t,x)'^ . (4.23) 

Consequently, (j4.10p reduces to the following system of two equations 

dtp^ + d,p+ = -[^0 + bil - uf]p+ + [^0 + b{l + uf]p- , (4.24) 
dtp--d^p- = -[jo + b{l + uf]p- + [jo + b{l - uf]p+ . (4.25) 
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Equivalently, defining the mass density g and flux j by 

Q{t,x) ■.= p^{t,x) +p~{t,x) , j{t,x) ■.= p^{t,x) -p~{t,x) 

the system can be written in the hydrodynamic description as 



dtQ + dxj 
dtj + dxQ 

u{t, x) 



0, 



2[jo + b{l + u^)]j + 4bgu , 



X i 5o[£>](t) 



^^w(\x - z\)Q{t,z)dz ' 

0, xG5o[i?](t) 



with 



5olel(i) 



w{\x — z\)Q{t^ z) dz = 



(4.26) 

(4.27) 
(4.28) 

(4.29) 



Remark 2 Due to the assumption [A2']\ we have p'^{t,x) = p {t,x) = on So[p'^,p ]{t). Therefore, 
both q^'^{x; x,u) in ( j^.i^P and r^'^{x; x,v) in l\4-15\) are equal to zero by definition. Consequently, 

q^{x;x,v)= lim q^'^{x; x,v) = , and r^{x;x,v)= lim r^'^{x; x,v) = , 
N^oo TV— >-oo 



and formulas l \4-S3\) remain valid irrespective of the particular choice of the value ofu{t,x). This justifies 
our extension of the definition of u{t,x) by setting it equal to zero for x S SqIp'^ , p~]{t) . 



5 Existence of solutions to the kinetic system 

The main goal of this Section is to prove the following Theorem: 



Theorem 1 Let > 0, b > and w satisfy the assumptions [Al]\ and [A2']\ Then, for every T > 
and every nonnegative initial datum {pq,Pq) G L°°(il) x L°°(f]) there exists a nonnegative solution to 
the kinetic formulation l \4-24\) - ^4-^5\ l L°°{[0,T] x 0). This also establishes solutions Q = p^ 
j = p^ — p~ of the hydrodynamic formulation l\4-27\ i - l\4-29\ i with the corresponding initial condition. 

Proof of Theorem [1} The proof is carried out in three steps and is only sketched here, omitting details 
where the techniques are standard. For notational convenience, we will work both with the kinetic and 
hydrodynamic representation of the system and treat {p'^,p~) and {Q,j) as synonyms, related by ()4.26p . 

Step 1. First, we consider a linearized version of (j4.24p ~ (j4.25p . where we solve for p'^ and p~ given a 
prescribed u G L°°([0,T] x il) with \u\ < 1. This constitutes a strictly hyperbolic system with unique 
mild, nonnegative solution in C{[0,T]; L°°{Q,)), constructed by a standard fixed point iteration (see, for 
instance, [10], Section 7.3). The solution is given by the Duhamel formula 



p'^{t,x) 
p~{t,x) 



p+{0,x-t)+ [ Q+[p+,p~]{s,x + {s-t))ds. 
Jo 

p~{0,x + t)+ / Q'[p^,p'~]{s,x-{s-t))ds. 
Jo 



(5.1) 
(5.2) 
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with Q^\p^ ,p~] = T[lo + b{l — lb [70 + 6(1 + Moreover, for any fixed T > 0, we have apriori 

boundedness of p'^ and p~ in L°°{[0,T] x O), depending only on the initial condition. Indeed, denoting 
h(t) := sup^gQp+(i, x) + sup^^Q p~{t,x), and remembering the uniform boundedness \u\ < 1, we have 

h{t) < h{0) + (70 + 46) / h{s) ds , 

Jo 

and the apriori boundedness follows from an application of the Gronwall lemma on the time interval [0, T]. 
Step 2. We consider a regularized version of (|4.24p - ()4.25p where u is substituted by n^, defined by 
Jix) 



Ue{x) :-- 



e + R{x) 



with J{x) := / w{\x — z\)j{z) dz , R{x) := / w{\x — z\)q{z) dz . (5.3) 
Jn Jn 



For any fixed e > 0, a solution is found by the Schauder fixed point iteration on the mean velocity Us [10] . 
The compactness of the corresponding Schauder operator is provided by the Arzela-Ascoli theorem. In- 
deed, let us take a sequence it"" with ||'w"||j;,oo([o T]xf7) — ^ P~^'^ P~'^ be the corresponding 
mild solutions of the kinetic system, given by the Duhamel formula (j5.ip - ()5.2p with in place of u, 
and let = p"*"'" + and j" = p^'" — As explained in Step 1, for any fixed T > we have 
IIp^'^II^cxiqo T]xn)i IIp^'"IIl°°{[o T]xn) bounded uniformly with respect to n. Defining the function 

oj{x) := I \w{z) — w{\z — x\)\dz for x G , 

one has lima;_^o ^^(2;) = (continuity of translation [21]). By Holder inequality, 

|J"(x)-J"(y)|<||r||^^(^)a;(x-y), 
with J'^{x) := j^w{\x — z|)j"(z)dz. Moreover, we have the uniform boundedness 



iL°°(n) 



\W\\t.i 



and analogous estimates hold for R^{x) := J^w{\x — z\)Q{z)dz. Consequently, we have 

< ^11/ 

e 



This uniform equicontinuity together with the uniform boundedness < 1 allows us to apply the 
Arzela-Ascoli Lemma and obtain the compactness of the Schauder operator. Therefore, for every fixed 
e > 0, we have a nonnegative solution (pf,pj) of the regularized system (j4.24p . (j4.25p . (j5.3p . uniformly 
bounded (with respect to e) in L°°([0,T] x $7). 

Step 3. Finally, we pass to the limit e — > 0. Due to the uniform boundedness, a subsequence of pf p^ 
weakly* in L°°([0, T]; L°°(ri)); we need to show that the limit of the nonlinear terms QgU^ and u^je is gu 
and, resp., u^j, with u given by (j4.22p . or, equivalently, (j4.29p . The limit passage in the distributional 
formulation of the term g^u,; with a test function (p G C^([0, T) x 17) is performed as follows: 



roo r rco r roo r 

/ / {geUe - gu)Lpdxdt < / {ge - g)uLp dx dt + 

Jo Jn Jo Jn Jo Jn 



ge{ue — u)ipdxdt 
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The first term vanishes in the hmit e ^ due to the weak* convergence of Qs towards g with mp a valid 
test function. Concerning the second term, we will show that Qeiue — u)ipdx tends to zero for almost 
all t € [0, T] and conclude the convergence of the time integral by the Lebesgue dominated convergence 
theorem. Let us fix t G [0, T] and define the set Ss by 

Ss := {x G n; R{x) < 5} , 

with 5 > and R defined in ()5.3p . Then, we have 



/ Qeiue - u)(pdx < 2 / 
JSs JS 



< 2 I ^e|v:'|dx > 2 / Q\ip\ 



dx 



= 2 / Q\ip\ dx , 

Jss\So <5->o 

where the second line is due to ^? = on 5o and because meas(55 \So) tends to zero as 5 — )■ 0. Next, for 
X e ^\Ss, 



R{x)Re{x) 



R{x) 



< -{\Re{x)-R{x)\ + \Je{x)-J{x)\), 

where Je(x) = J^w{\x — z\)ji;{z)dz and R£{x) = f^w{\x — z\)qs{z) dz. Therefore, due to the uniform 
convergence of i?e and to R and, resp., J onil. (implied by the uniform equicontinuity and boundedness 
of the families {Re}e>o and {Je}e>o), we have 



/ 



QsiUs - u)ipdx 



< I supf^\5jn, - m 



Qe\^\ dx 



~ ^ v'-^^ ~ 11"'^ ~ ll'^'lli'Mf^) 



We conclude by passing 5^0. The limit passage in the term u^je is performed similarly (note that 

\je\ < Qe)- 



It is worth noting that the assumption | [Al] | of Theorem[T]can be relaxed. In fact, we posed the requirement 
[Al] of boundedness of w on VL in order to establish the definition ()2.2p of uf" in the discrete model and 
pass to the limit N ^ oo. However, at the level of the kinetic or hydrodynamic description, we may 
relax this to w (z L^(r2). Moreover, formally it is possible to consider even singular weights, in particular 
w = (5o, which leads to u = j/ q and removes the nonlocality. In fact, one can see the choice w = 5q as the 
limiting case when the interaction radius shrinks to zero: for almost all x G such that q{x) ^ 0, one has 



lim 



Lx[o,<7](|a;-^;|)j(2)d2; _ i{x) 



"^^o !n X[o,<7] {\x - z\)q{z) dz q{x) ' 

where X[o,(t] is the characteristic function of the interval [0,cj]. One can interpret this as a model where 
only pointwise local observations of the system are possible. 
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By a slight modification of the proof of Theorem [T] it is possible to show that given a sequence of 
weights Wn converging strongly in L^(il) to w, the solutions {pn,jn) corresponding to the weights Wn 
converge weakly* in L°°{[0,T]; L°°{il.)) to the solution corresponding to the weight w. However, we need 
the condition [A2'] to be satisfied uniformly; consequently, the question whether and how the solution 
corresponding to w = Sq can be obtained as a limit of solutions with Wa = X[o,a] as it — )• remains an 
interesting open problem. An even more interesting question is what is the limit of the discrete model as 
— 7> oo if the interaction radius is shrinking as some power of 1 /N. This question is studied in Section [7] 
below. 

6 Long time behaviour 

In this Section we provide several conjectures about the long time behaviour of the kinetic system (|4.22p . 
(j4.24p . (j4.25p or, equivalently, the hydrodynamic system (j4.27p - (l4.29p . To get some insight into the long 
time dynamics, we start with a numerical example. We solve the kinetic system using standard semi- 
implicit finite difference method with upwinding. The initial condition is Pq = 2.2 on [0.125,0.375] and 
zero otherwise, pg = 1.8 on [0.625,0.875] and zero otherwise. In Figure [6] we show the results for the 
choice of parameters 6 = 1 and 70 = 0.3 and the weight function w = X[o,o.2]- 

We conjecture that with regular weights w satisfying [Al] and [A2'], the solutions {p,j) to (I4.27p - ()4.29p 
converge to the constants /O = 1 and j = js with some js G K, exponentially fast as t — )• 00. Moreover, 
in the large noise case 70 > b, we hypothesize that js = 0. Unfortunately, we are only able to provide an 
analytical proof in the rather special case w = 1 and 70 > b: 

Lemma 3 Assuming w = 1, we have u{t,x) = u{t), where u(t) satisfies the ordinary differential equation 

n = -2(7o + 6(n2-l))u, (6.1) 
subject to the initial condition u{0) = J^j{0,x)dx. Moreover, 

(i) if 70 < b, then Y\m.t^^u{t) = sign(u(0))A/l - 706-^, 

(a) if Jo > b, then limt^oo u{t) = and 

< h(0)|e-2(^«-'')*. (6.2) 
Moreover, j converges to zero exponentially fast in the L"^ -sense: 

I j2(t,x)da; < ce-^To* 
Jn 

for a suitable constant c. 

Proof: Integrating (j4.28p with w = 1 with respect to x € we obtain (j6.ip . This has the stationary 
state u = 0, which is stable if and only if 70 > b. Moreover, if 70 < b, two additional stable stationary 
states u = iby^l — 7o6~^ exist. This establishes the first statement. Further, we have 

^ Qn^) = -2(70 + biu' - 1)W < -2(70 - bW , 
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Figure 3: Numerical results with w = X[o,o.2]j b = 1 and 70 = 0.3: and p~ converge to constant states, 
the mass density g converges to 1, the flux j (full line) and averaged velocity u (dashed line) converge to 
a constant. 
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and an application of the Gronwall lemma gives (j6.2p . 

To prove the convergence of j to zero, we consider the identity 



1^ 
2 dt 



/ (e*^ + j^) dx = -2(70 + 6(1 + u^)) / fdx + ibu / gjdx. (6.3) 
Jn Jn Jn 

An application of the Cauchy-Schwartz inequality and the decay rate (j6.2p yield 

Abu I gjdx<2b [ fdx + 2bu^ [ g"^ dx < 2b [ dx + 26|w(0)pe-^(^o-^)* [ g"^ dx . 
Jn Jn Jn Jn Jn 

Then, an integration of ()6.3p in time leads to 

[ ^2(r,x)dx < co + 46|^/(0)pe-^(^«~^)* /" [ g'^it,x)dxdt , 
Jn Jo Jn 

with Co := gQ{x) + jQ{x) dx. Consequently, by the Gronwall lemma, g^ dx is bounded uniformly in 
time by a constant ci if 70 > b. Inserting this information into (|6.3p gives 

[ f{T,x)dx < co-470/" [ f{t,x)dxdt + ibci\u{0)\'^ [ 6-^(^0-'')* dt 
Jn Jo Jn Jo 

< c-470 /" [ f{t,x)dxdt 
Jo Jn 

and an application of the Gronwall lemma yields the second statement. 



6.1 The case w = 60 

In this subsection we briefly discuss the long time behaviour in the singular case w = 5o. Again, we start 
with a numerical example where we solve the kinetic system with the parameters b = 1 and 70 = 0.3. The 
initial condition is chosen as before, see Figure [6] (left panels). In Figure [6TT] we present the time evolution 
of p~^, p~ p, j and u. Based on the numerical observations, we conjecture that, for the small noise case 
7o < b, the long time dynamics are given by the travelling wave profiles {pf,pj), satisfying 

dtpt + dxPt = , 
dtP'^ - dxPj = , 



and ^1 G {uc,, —Ug} with Us = \/l — 706 Then, it is a matter of a simple consideration to deduce 

Ps ~^Ps 

that one of the functions, say pj , has to be a global constant, while the other one, pf, is a piecewise 
constant assuming only two values {p^i,pt2}^ satisfying the relations 



P~^ / 1 w \ ^ 
Pt,pt2 = iPsf and = 



These relations chracterize the dynamic equilibria between the densities of individuals marching to the 
left and to the right, in dependence on the parameter values 70 and b. 
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Figure 4: Numerical results with the singular weight w = 5q, b = 1 and 70 = 0.3: p~ converges to 
a constant, while becomes a piecewise constant travelling wave profile and u (dashed line) jumps 
between the values ±Us with Ug = \/l — job~^ = 0.8367. 
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For the large noise case (70 > b), we prove that the asymptotic state is again i = and p = 1. In 
fact, j converges to in the L^-sense exponentially fast as t — > 00. This follows from an application of 
the Gronwall lemma to the estimate 

< -2(70-6) j fdx. 

7 The effect of shrinking interaction radius 

In biological applications, the weight function w is often chosen as the characteristic function of the 
interval [0, a], where a > is the interaction radius. It is important to study the dependence of collective 
behaviour on the size of o". In this section, we will consider a theoretical limit o" — )■ 0. Clearly, letting a ^ 
with fixed N leads to a trivial model without any interactions (it corresponds to the choice b = 0). On 
the other hand, letting first N — ?• 00, followed by o" ^ 0, we obtain the model with w = 6q, as mentioned 
in Section [5l Consequently, the limits o" ^ and — )• 00 do not commute. From the modelling point 
of view, it is natural to study the limit N ^ 00 with the interaction radius shrinking as A^""^ for some 
a > 0. Indeed, as the space is getting more crowded, the visibility is reduced and every individual can 
only take into account its closest neighbours. Actually, taking w = X[o,N-'^] letting — )■ 00, we will 
show that a significant limit is obtained with a = 1, leading to a new kinetic model. The results are 
summarized in the following Theorem. 



Theorem 2 Let w = X[o,n~"]- Then, the formal limit as N ^ 00 of the BBGKY hierarchy d^.iO ) is 
(1) For a > 1: 



(2) ForO<a< 1: 



(3) For a = l: 



d''tQ + 2{jo + b)dtQ = dlQ. (7.1) 



dtQ + d.^3 = 0, (7.2) 
dtj + d^Q = -2(^7o + 6(^l + ^))i + 46j. (7.3) 



dtQ + d^j = 0, (7.4) 
dtQ-d.,j = -2(7o + 6(l + r?))j + 46j(l-exp(-^)), (7.5) 



with 



7] = — {I- exp(-^)) + (1-^1 exjp{-g) [Et{Q) - 7 - ln{Q)] , (7.6) 

where Ei(z) = fl^ exp(— s) ds is the so-called exponential integral function and 7 = 0.577 ... is 
Euler's constant. 
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We see that for < a < 1, we obtained the kinetic model with u = j/p (i.e., the same as if one would 
first pass to — t- oo and then o" — )> 0). If q > 1, we get the model with no interactions (i.e., as if 
one would first pass to (j — )> and then — )> oo) which is described by the telegraph equation for g 
|19j . In the significant limit with a = 1 we obtained a new model, which is in fact the hydrodynamic 
model (|4.27p ~ (|4.28p with the function ^(1 — exp(— £>)) in place of u and with rj, given by (j7.6p . in place of 



Proof: All we need to do is to recalculate the limits A'" — )• cx) in the expressions (j4.1ip and (I4.12P with 
w = X[o,N-<^]- Let us fix z € and assume that it is a Lebesgue point of g and j, i.e., 

1 r 

lim - / g{z - y) dy = g{z) , (7.7) 



o— 5-0 a jQ 

and similarly for j. Moreover, assume that g{z) > 0. In the same way as in Lemma [H we calculate 

„ V^7V-l /I IN N-1 

^ ^ f Z^m=lX[0,N—]{\z - Xm\)Vm -pr 

Qn := 2^ / — ^7v3l 7 — [[p{xi,v^)dJc 



1 



Jn^-^ l + EZ=iX[o,N-'']i\z-Xm\) 
where we denoted 

In ■■= / X[o,N-'']iy)Q{z - y)dy , Kn := / x[o.Af-"] (y)i(^ - y) dy . 

Obviously, the sum X^^Zj X[o,Af-'»](k ~ s;^!) only takes the values A; = 0, . . . , A^ — 1, and 



N-2-k 



Therefore, 



/ 



I n \ / 



f^^-^ l + Em=IX[0,7V-](l 



■2' 



A,'=0 

1 



Due to (|7.7p . we have limTv-j-oo ^°'In = q{z) and limTv^oo N'^Kn = j{z). Therefore, 
hm Qr, = hm ^ (l - (1 - 1^)^^^) = ^4 (l - (1 - /tv)^"^) 



and 



hm (1 - (1 - /;v)^-^) = 



1 — exp(— ^(z)) for a = 1 , 

1 for < a < 1 , 

for a > 1 . 
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Using the same technique as above, one calculates that 

J-N \ -'at / 



with 



k=l ^ 



It is easily shown that Gn{In) vanishes in the limit N ^ oo whenever a ^ 1, while with a = 1 the law 
of rare events ( [H] ) gives 

lim Gn{In) = exp{-g{z)) V —- = exp{-Q{z)) [E[{q{z)) - 7 - ln(£>(z))] , 

Af— >oo ^ — ' Kl K 

k=l 

We assume that g and j be integrable functions, so that almost every z € is a Lebesgue point for them. 



Remark 3 Our analysis can he seen as a first step towards a so-called topological model, where the agent 
interactions are based on some connectivity graphs. For instance, one can consider the situation where 
every agent interacts only with its nearest neighbour. In this case, we have a different time dependent 
weight function Wi for every agent, namely Wi = X[o,o-,] with cjj = mium^j |xj — Xm\- Then, the passage to 
the corresponding continuum model is a completely open problem; however, let us observe that 



—y 

- 1 ^ 



mm \ Xi 



N — 1 ^ \xi — Xmr / p^-oo m 

Therefore, it would be interesting to consider the discrete system with Wi = X[o.o-i]; 



and study the limit — > 00 and p ^ 00 (possibly with p = N). Moreover, let us observe that "on 
average", \xi — Xm\ ~ . Therefore, 



. -i/p 

N - I ^ \xi - a 



so we are in the situation of the significant limit a = 1 described above, and we might believe that the new 
kinetic model obtained in this limit can have some connection with the limit N ^ 00 and p ^ 00 of ( |7.<§D . 
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8 Discussion 



We introduced an individual based model with velocity jumps aimed at explaining the experimentally 
observed collective motion of locusts marching in a ring shaped arena [2]. The frequency of individual 
velocity jumps increases with a local or global loss of group alignment. We showed that our model has the 
same predictive power as the model of Czirok and Vicsek, in particular, it exhibits the rapid transition to 
highly aligned collective motion as the size of the group grows and the switching of the group direction, 
with frequency rapidly decreasing with increasing group size. Moreover, in the limit — > oo we obtained 
a system of two kinetic equations. We proved existence of its solutions and a partial result about the 
long time behaviour. Finally, we studied the effect of shrinking the interaction radius a in the discrete 
model as the number of individuals, A^, tends to infinity. We showed that in the significant limit where 
a shrinks as one obtains a new kinetic model. 

Kinetic approach has previously been used in the literature to understand collective dynamics of 
individual based models. Carrillo et al [3] found the double milling phenomena in the kinetic formulation 
of the model of self propeled particles with three zones of interactions. The kinetic description of the 
Cucker-Smale model was introduced in [M], which can also be be derived from the Boltzmann-type 
equation, see H], or Povzner-type equation, |12j . For the survey of the most recent results see the 
review [5]. 

Several interesting questions remain open, offering space for future investigations. For example, the 
kinetic system (j4.24p - (j4.25p deserves a better analysis, in particular, uniqueness of solutions and more 
complete investigation of the long time behaviour. It would also be interesting to know if and how the 
solutions corresponding to w = Sq can be derived as a limit of solutions corresponding to w = X[o,a] 
(T — 7> 0. Another interesting direction of future research was formulated in Remark [3l 
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From individual to collective behaviour of coupled velocity 
jump processes: a locust example 



Radek Erban^ Jan Haskovec^ 



Abstract. A class of stochastic individual-based models, written in terms of coupled velocity jump 
processes, is presented and analysed. This modelling approach incorporates recent experimental findings 
on behaviour of locusts. It exhibits nontrivial dynamics with a "phase change" behaviour and recovers 
the observed group directional switching. Estimates of the expected switching times, in terms of number 
of individuals and values of the model coefficients, are obtained using the corresponding Fokker-Planck 
equation. In the limit of large populations, a system of two kinetic equations with nonlocal and nonlinear 
right hand side is derived and analyzed. The existence of its solutions is proven and the system's long-time 
behaviour is investigated. 

Key words: Collective behaviour. Stochastic individual-based model. Density-dependent directional 
switching, Kinetic equation. 

1 Introduction 

Individual-based behaviour in biology can be often modelled as a velocity jump process [20]. Here, the 
velocity of an individual is subject to sudden changes ("jumps") at random instants. If the velocity changes 
are completely random, this process simply leads to diffusive spreading of individuals in an appropriate 
limit [18]. The situation becomes more complicated whenever the velocity changes are biased according to 
an individual's environment. A classical example is bacterial chemotaxis [8]. Individual bacteria change 
their frequency of velocity changes according to their environment. If they swim in a favourable direction 
(e.g. towards a nutrient source), they are less likely to change their direction. On the other hand, they 
are more likely to turn if they are heading away from a foodstuff [9] . 

In this paper, we modify the velocity jump methodology to model the behaviour of locusts. Our model 
is motivated by the recent experiments of Buhl et al [2]. They studied an experimental setting, in which 
locust nymphs marched in a ring-shaped arena. The collective behaviour depended strongly on locust 
density. At low densities, there was a low incidence of alignment among individuals. Intermediate densities 
were characterized by long periods of collective motion in one direction along the arena interrupted by 
rapid changes of group direction. If the density of locusts was further increased, the group quickly adopted 
a common and persistent rotational direction. Yates et al [25] analysed experimental data of Buhl et al [2] 
and proposed that the frequency of random changes in the direction of an individual increases when the 
individual looses the alignment with the rest of the group. In this paper, we incorporate this observation 
into a stochastic individual-based model formulated as a velocity jump process. We show that this model, 
although phenomenologically very simple, has the same predictive power as other modelling approaches 
previously used in this area [7, 2]. In particular, it exhibits (i) a rapid transition from disordered movement 
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of individuals to highly aligned collective motion as the size of the group grows, and (ii) sudden and rapid 
switching of the group direction, with frequency decreasing as group size increases. 

The individual-based model is introduced in Section 2. The ring-shaped arena, used in Buhl's experi- 
ments [2], is modelled as one-dimensional interval with periodic boundary conditions. Locusts march with 
a constant speed and each individual switches its direction randomly. The individual switching frequency 
increases in response to a loss of alignment. In Section 3, the corresponding Fokker-Planck equation is 
derived for the system with global interactions and possible types of qualitative behaviour of the system 
are classified. For the case of ordered group motion, where two distinct metastable states exist, an ap- 
proximate analytic formula for the mean switching time between these two states is derived. Then, in 
Section 4, the kinetic formulation of the model is obtained in the limit as the number of locusts tends to 
infinity. The existence of solutions of the kinetic model is shown in Section 5 and the long time behaviour 
is investigated in Section 6. We conclude with analysis of the dependence of collective behaviour on the 
size of the interaction radius of the individuals in Section 7. 



2 Individual based model 

We consider a group of N agents (locusts) with time-dependent positions Xi{t) and velocities i = 

1, . . . ,N. To mimic the ring-shaped arena set-up of [2], we assume that the agents move along a one- 
dimensional circle, which we identify with the interval = [0, 1) with periodic boundary conditions, and 
move either to the right or to the left with the same unit speed, i.e. 

Xi{t) en, Vi£ {-1, 1} and -^{t) = Vi{t). (2.1) 
We define the local average velocity of the ensemble, seen by the i-th. agent, as 

,loc _ YlZ=iw{\xi -Xm\)u 

YlZ=lW{\xi-Xm\) 

where w is a weight function defined on Q with the properties: 
[Al] w is bounded and nonnegative on 0, 
[A2] w{0) > 0. 

For example, w(s) = X[o,o-]('5)> where X[o,cr] is the characteristic function of the interval [0, a] and o" > 
is a interaction radius, satisfies conditions [Al] and [A2]. This is a common choice of w in biological 
applications [7, 25]. It is worth noting that, due to the assumption [A2], the definition (2.2) always makes 
sense and u{°'^ G [— !]• 

The agents switch their velocities to the opposite direction (i.e., from Vi = 1 to Vi = —1 and vice 
versa) based on independent Poisson processes with the rates 

li = 7o + bC{vi-u'f^), i = l,...,N, 

where 70 > and 6 > are fixed parameters and the "response to disalignment" function ^ : [—2,2] — )> 
[0, 00) is assumed to be convex, differentiable and symmetric with respect to the origin. Taking the Taylor 
expansion of ^(s) around s = 0, we obtain 

C(s) = ao + "25^ + C(s^) • 
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N=5 N=7 N=12 




Figure 1: An example of the transition to ordered motion as N grows. We use (2.1)-(2.4) with 6=1, 
7o = 0.2 and w = X[o,o.2] • Shown are the normalized histograms of the group mean velocities u = jj "^iLi Vi 
recorded in 10^ time steps, with N = 5 (left panel), N = 7 (middle panel) and N = 12 individuals (right 
panel). The system does not prefer any particular state for N = 5. Two quasi-stable states of ordered 
collective motion are easily recognizable for N = 12. 

We can set = without loss of generality, because it can be absorbed in 79. Since the individuals 
switch their velocities less frequently when they are aligned ([25]), .^(s) has a global minimum at s = 0. 
This implies that 02 > 0. If 02 > 0, we can set a2 = 1 by choosing an appropriate time scale. Therefore, 
^(s) has the general form + 0{s^). For the rest of the paper, we choose the form ^{s) = s"^ for 
simplicity. Other choices are certainly possible, for example, in the limiting case 02 = the leading order 
approximation is given by a higher order term, which, however, complicates the analysis. However, it is 
worth noting that the derivation of the kinetic equation performed in Section 4 is possible, for example, 
also for ^(s) = \s\^, n > 3. 

With ^(s) = s^, the turning rate "from the right to the left", J^^^, and the rate for the opposite 
turn, 'y^^^^ are given by 

^R^L ^ + 5 (1 _ ^«oc^2 ^j^g switch from t;^ = 1 to 1;^ = -1 , (2.3) 
7/-^^ = 7q + 6 (1 + u^°c^2 ^Yie switch from Vi = -1 to Vi = 1 . (2.4) 

This velocity jump process describes the tendency of the individuals to align their velocities to the 
average velocity of their neighbors. Despite its relative simplicity, the model provides similar predictions 
as the Vicsek and Czirok model [7] and its modification [25] and is in qualitative agreement with the 
experimental observations made in [2]: the transition to ordered motion as grows (Figure 1) and the 
density-dependent switching behaviour between the ordered states (Figure 2, bottom). 
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Figure 2: An example of the behaviour of the model (2.1)~(2.4) with global interactions (w = 1): the 
large noise case (top) with N = 20, 6=1 and 70 = 1.3 and small noise case (bottom) with N = 20, 
6 = 1 and 70 = 0.3. Left are the histograms of the group mean velocities u recorded in 10^ time steps, 
compared to the plot of the (properly scaled) stationary solution ps of the corresponding Fokker-Planck 
equation (solid line). Right are the plots of the temporal evolution of the group mean velocity u during 
5 X 10^ timesteps. In the small noise case, one can clearly distinguish the two quasi- stationary states and 
observe the switching between them. 



4 



3 Analysis of the individual based model with global interactions 

In this section, we simplify the individual-based model by assuming w = 1 in (2.2), i.e. u[°'^ = u for all 
i = 1, . . . , N, where 

1 ^ , , 2r(t) -N , , 

4 = 1 

where r(t) is the number of individuals which are going to the right at time t. Using this simplification, 
we will obtain an explicit formula for the mean switching time between the ordered states. However, 
for the derivation of the kinetic decription and its analysis (Section 4), we will allow general weights w, 
imposing only the assumption [Al] and a slightly reinforced version of [A2]. 

Let p{r,t) be the probability that r individuals move to the right (i.e., with velocity 1) at time t>0. 
It satisfies the master equation 

^p(r,t) = (r + l)7«^^p(r + l,t)-r7^^V^t) 

+ {N -r+l)'y^~^^p{r -l,t) - (N -r)j^-^^p{r,t) , (3.2) 

where and 7^^^ are given by (2.3) and (2.4), respectively. The subscript i in (2.3)-(2.4) is dropped 

in (3.2) because all individuals have the same turning rates. Using the system size expansion [24] and the 
definition (3.1) of the average velocity u{t), we obtain the following Fokker-Planck equation 



dp{u, t)_d_ _ _ ^^^^ ^ 1^ ^ _ ^^^^ 



dt 



(3.3) 



where p{u, t) is the probability distribution function of the average velocity (3.1) at time t. The stationary 
solution ps of (3.3) is 

Ps{u) = Cexp[-^Niu)] (3.4) 
where C is the normalization constant and the potential <&Ar is given by 

^n{u) = -^u2 + (1 - ^) In (70 + 6(1 - n2)) . (3.5) 

The comparision of the stationary probability distribution function ps with the results obtained by long- 
time simulation of the stochastic individual-based model is shown in Figure 2. Differentiating (3.5) twice, 
we obtain 

Consequently, we distinguish the following two cases: 
7q 2 

(1) Large noise: If > 1 + — , then has the global minimum at u = 0. 

(2) Small noise: If ^ ^ "f" ' ^^^^ ^ local maximum at u = 0. The only local and global 
minima are at it n, where 



2 70 
IH -. 

N b 



5 



It is worth noting that > 1 if is smah, namely for < 26/79. This is a consequence of approximations 
made during the derivation of the Fokker-Planck equation (3.3). 

In the large noise case, the system prefers the disordered state u = 0, while, in the small noise case, 
the system has two preferred ordered states ±u<j where it spends most of the time and eventually switches 
between them (see Figure 2, bottom panel). 

Using Kramers theory [15, 13], the mean switching time tn between the states —Us and Ug can be 
approximated as 

f ^ exp($jv(0) - ^jv(-n,)) 

Tm[—Us i-> Us) ~ r = — , 3.7) 

^ ' 70 + 6 V-<I>'^(0)«I>'^(-n,) ^ ' 

which has an exponential asymptotic with respect to large N given by 



b--foy 2(70 + b) 



&-70 _ To 270 



70 + 6 



This is in agreement with the experimental observations, [2], as well as with the modified Czirok-Vicsek 
model of [25], where the mean switching time is as well exponential in A^. Finally, it is interesting to note 
that with the transform b = 796 in (3.7), one has 

r(iV,7o,6) =7or(A^,6), 

i.e., if 6/70 is kept fixed, the mean switching time scales linearly with 70. 

In the numerical experiment shown in Figure 2 bottom (small noise case with 70 = 0.3, 6 = 1 and 
= 20 agents), we have two metastable states located approximately at Us = ±0.894. The estimate 

mean turning time given by formula (3.7) is tn{—Us ^ Us) = 61.1. Performing 10^ time steps of length 

10~^, the observed mean switching time (defined as a mean transition time between the states v = —0.8 

and V = 0.8 or vice versa) was 58.8, showing a very good agreement. 



4 Kinetic description 

In this section we derive the kinetic description of the system of A^ interacting agents and formally pass 
to the limit A^ — > 00 to obtain the corresponding kinetic equation. For this, we have to accept a slight 
reinforcement of the assumption [A2] on namely, 

[A2'] > on the interval [0, r) for some r > 0. 

The state of the system of A^ agents at time t > is described by the probability density function 
p^(t, . . . ,X]\i,V]\i) of finding the i-th agent in position Xi ^ Q with velocity Vi G {—1, 1}, for i = 
1, . . . ,N. When convenient, we will use the abbreviation = p^{t, x, v), with x = (xi, . . . ,xn) G 
and V = {vi, . . . ,vn) G {±1}^. The probability density p^(t,x, v) evolves according to 

|^'^ + E-^^^'^ = -^[^'^] + ^b'^]' (4-1) 

1=1 * 

where i2[p''^] (resp. ^[p^]) is the loss (resp. gain) term corresponding to the velocity jumps. Using 
(2.3)-(2.4), the rate of the switch Vi —Vi is for z G {1, . . . , A^} given by 

7i(x, v) = 70 + 6 (tij - u[xi;x,v])'^ , 
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where u[2;j;x, v] = v}°^ is defined by (2.2), i.e. u[xj;x, v] is the local average velocity seen by an agent 
located at Xi € $7, based on the system configuration [x, v]. Consequently, the loss term is given by 



TV 



/:[p^](f,x,v) = J]7,(x,v)p^(t,x,v). (4.2) 

i=l 

We define the operator Mi : {±1}*= {±1}'' for A; € {1, ... , N} and z G {1, . . . , k} by 

Mj(v) = {vi, . . . ,Vi-i, -Vi,Vi-i, . . .,Vk), (4.3) 

i.e. Mj(v) denotes the velocity vector created from v by changing the sign of its i-th component. Then 
the gain term Q[p^]{t,'x.,v) is given by 

TV 

g[p^](t,x,v) = J]7,(x,M,(v))j>^(t,x,M,(v)), (4.4) 

i=l 

where 

7j(x,Mj(v)) = 70 + & (- t-i - u[xi;x,Mi{v)])'^ . 
Consequently, the right hand side of (4.1) is 

(g-/:)[p^](x,v) = 7o(^J^p^(x,Af,(v))-iVp^(x,v)^ (4.5) 

TV 
1=1 

where we dropped the dependance on time t to simplify the notation. Finally, we postulate the so-called 
indistinguishability-of-particles: We only consider solutions p^ that are indifferent to permutations of their 
(xj, Uj)-arguments. Such solutions are admissible, since the equation (4.1) with the collision operator (4.5) 
is as well indifferent with respect to interchange of the (xj, 'L'j)-pairs. 

4.1 Derivation of the BBGKY hierarchy 

To derive an analogue of what is called the BBGKY hierarchy in the classical kinetic theory of gases (see, 
for instance, [6]), we define, for k = 1, . . . , N , the fc-agent marginals 

pN,k^^k^^k^^ j2 I p^(x^x,v^v)dx, x'^ g o^ g {±l}^ (4.6) 



VG{±1}^ 

In what follows, the coordinates of x'^ and x will be denoted as 

x'^ = {x\, X2, . . . , X^) , X = (xi, X2, . . . , Xj^-k) , 
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ve{±i} 



that is, X = (x^,x) = (x^;', x^, • • • , a;^, xi, X2, . . . ,^Ar_fc). The same notational convention is used for 
velocities, i.e., v = (v'^', v) = {vi,V2, ■ ■ ■ ,v^,vi,V2, ■ ■ ■ ,VN~k)- Note that, due to the indistinguishabihty- 
of-particles, the marginals are well defined and indifferent with respect to permutations of the pairs of 
arguments {xi,Vi), i = 1, . . . , k. Integrating (4.1), we obtain 

|p^''=(x^v'=)+^.,Ap^'.'=(x^v'=)= L /a-^)[p^](x^x,v^v)dx. (4.7) 

Substituting (4.5) into the right-hand side of (4.7), we get 

/^_^(a-/:)b^](x^x,v^v)dx 

k 

= 70 5](p^•'=(x^M,(v'=))-p^''=(x^v'^)) 

i=l 

+ ^ E ^E{(-^*'-^[^-^^^'^*(^')'^])V(x^x,M,(v^),v) 

- (ff - u[xi;x'^,x, v'^, v]) V^(x'^,x, v'', v)| 

k 

= (70+6) ^(p^''=(x^M,(v'=)) -p^''=(x^v^)) 
1=1 

+ 26 E / (n[x,;x^x,M,(v^•),v]p^(x^x,M,(v^),v) 

V6{±1}~-''^' *=1 

+ u[xj;x'^,x, v'^^ v]p^(x''\x, v'^, v)^ dx 
+ ^ E ^E((^[^-^''^'^^(^')'^])'p'^(^''X,M,(v'=),v) 

- (n[x,;x^x,v^v])'p^(x^x,v^v)) dx. (4.8) 



Since we are interested in the limit N oo with k fixed, we can rewrite the definition (2.2) as follows 



n|Xi;x ,x, V ,vj — ^ j^—j^ — 

Em=l ^(l^^m + Ylm=l ^(l^m " Xi\) 

-- n[xi;x,v] + 0( 41 , (4.9) 



N 



where we define 

yN-k 



n[z;x,v]= ^-^"^'"--"l^^- for.EO. 

22ra=lWi\^rn- Z\) 
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Substituting (4.9) in (4.8) and (4.7), we obtain the BBGKY hierarchy 



i=l ' 1=1 



+ 2b^v^ (q^'''{xi; x^Mi(v'=)j +q^'''(^Xi; x^v'=jj (4.10) 

1=1 

+ (^^; x^M,(v'=)) - r^'^(x.; x^ v'^)) + O (A) , 



where 



ve{±i} 



/ 'u[z;x,v]p^(x^x,v^v)dx. (4.11) 



and 



z; x'^jV^j = / 'u[z;x, v]^p^(x^',x, v'', v) dx. (4.12) 



^e{±i} 



4.2 Passage to the limit N ^ oo 

The usual procedure of deriving the mean field equation is to write the BBGKY hierarchy (4.10) in terms 
oip^'^ and pass to the limit ^ oo to obtain the so-called Boltzmann hierarchy for p'^ := limAr^ooP^'*^ 
[23]. Then, one shows that the Boltzmann hierarchy admits solutions generated by the molecular chaos 
ansatz (see below). In our case, however, this strategy cannot be pursued; although we could derive 
uniform estimates allowing us to pass to the limit — )• cx) in the BBGKY hierarchy, we are not able to 
express the limiting 

^N,k terms of p^. Consequently, 
we have no clue what the correct molecular chaos ansatz for and r'^ should be. 

Instead, we assume the propagation of chaos already at the level of the BBGKY hierarchy, before 
passing to the limit N ^ oo: we assume that, for large N, is well approximated by the product of the 
limiting one-particle marginals p '. — liniTv^QQ p^'\ i.e. 

N 

p^{t,x,v) ^Ylp{t,Xi,Vi) for ah t > 0, x e Jl^, V e {±1}^. (4.13) 

i=l 

This corresponds to vanishing statistical dependence (correlations) between the agents as ^ oo and 
is the usual phenomenon observed in systems of interacting particles, see for instance [6] in the context 
of classical kinetic theory or [16, 17] in the context of biological systems. Moreover, if one interprets 
qN,k (^j-ggp^ ^N,k^ ^l^g £j,g|. (^pggp^ second) order moment of p'^ with respect to tt[z;x, v], then one can 
understand (4.13) as the moment closure assumption for the non-closed system of moments generated 
by (4.1). 

The essential point is that now we may insert (4.13) into (4.11) and (4.12) to obtain explicit expressions 
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for and in terms of p\ 

g^(z;x^v'=)= lim q^^^iz'^^^ 

N-^oo \ 



lim q'^'^iz; x'',v'' = (4.14) 



and 



"^(z;x^v'=)= lim r^'^fz;x^v^' 

np(x^.f)') X lim 5: / f^fe%^%T np(-^,^odx. 



r'^(z;x^v'') = lim r^"''^ z; x^vM = (4.15) 



\±x'^ jv^oo ^ QN-k \ V^^"" infix - zl 

Study of the limit iV ^ oo in (4.14) and (4.15) 

We start by setting k = 1, which is the case considered in Section 4.3. To simplify the notation, we drop 
the bars over x and v, and, without loss of generality, choose z = 0. First, we explore the symmetry of 
the expression (4.14) as follows: 

= yZ Tl^T — r^^—r^iPi^niA)-p{xm,-l)] V TT p(2;i,Vi)dx 



' Jf^iv-i (iV - l)57v(x) , , 

/ '^nc / J i^m) n £'(xi) dx 

Jn^-i [N - l)5Ar(x) .-^-L 

m=l / \ / j^m 



with the notation 



£»(x) := p(x, 1) +p(x, -1) , 
j(x) := p(x, 1) -p(x, -1) , 



and 



N-l 

Due to the normalization of p^, we have f^Q{x)dx = 1. Consequently, in what follows we denote by 
Pg(t) the time dependent probability measure corresponding to the probability density Q{t). Since w is 
bounded and nonnegative by assumption [Al], it is integrable with respect to Pg and we may define 

/ := / w{x)dPg{x) > 0. (4.17) 
Jn 
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The forthcoming analysis will be performed given the assumption / > 0; the case / = will be discussed 
in Remark 2. 

Lemma 1 Let Pg he a probability measure on $7 with density g and j £ L^{Q) such that \j\ < p. Let 
w : Q —?■ [0, oo) with w € L°°{Q) be such that the integral I defined by (4.17) is positive. Define 



Qn '■- 



w{xi)j{x 



S'Ar(x) 



1=2 



Then 



lim Qn = j[ w{y)j{y)dy. 



Proof: We can treat w{y) as a random variable with respect to the probability measure Pg{y)- The 
essential tool of the proof is the law of large numbers, which states that Sn{^) converges to / in measure, 
in the sense that for each e > 0, 



hm Pf-i ({x G 0^-1; |5;v(x) - I| > e}) = , 



N-^oo 



(4.18) 



where P^ denotes the {N — l)-fold tensor product of the probability measures Pg. Moreover, the 
existence of the m-th order moment of w with respect to Pg, 

\w{y)r dPg{y) < cx, , 

implies the rate of convergence (see [1]) 

lim {N - lY'-^P^-^ ({x G VL^-^- |5jv(x) - I| > e}) = . 

Let us denote by A]^{e) the set {x € fi^"^, |SAr(x) — /| < e} and by A'f^{e) its complement in Q 
Choosing < e < L/2 and x G 74jv(e), we have the estimate 



(4.19) 



w{xi) w{xi) 



Consequently, 



5Ar(x 
w{xi) w{xi) 



2e , , 
< jiw{xi) . 



Aiv{e) 'S'Af(x) / 

On the other hand, the integral over ^^(e) is estimated with 
wixx) w{x\) 



dPf-l(x)<| / ^Xi)dPf"Hx) = y. 



Sn{x.) 



I 



dP^Hx)< / ^dPf-i(x)+ / 



w{xi) 
A'L{e) 'S'Ar(x) 



dP, 



N-l, 



X 



The first term on the right hand side converges to zero as A'^ — )• oo by (4.18) and the boundedness of w. 
Using (4.16), the second term is estimated by 



JA<i,(e) Sn{^) ^ 



\^)<iN-l)P'''\A%{e)) 
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and vanishes as ^ oo due to (4.19) with m = 2. Consequently, we have shown that the term 

N-l 



Qn -J [ w{y)j{y)dy 
^ Jn 



< 



w{xi) w{xi) 
QN-i \Sn{x) I 

wi^xi) w{xi) 



1=2 



I 



dPf-Hx) 



can be made arbitrarily small, for sufficiently large A^. 



The formula (4.15) can be analysed in a similar way as we did for (4.14). Namely, we exploit its symmetry 
as follows. 



E 



Jn^-r I (iV-l)5^(x) 



^2 N-~l 



Y[ pixi,Vi)dx 



VG{±1P 

(Af-1)2 ^ ^ J^N-l 



i=l 



w{Xi)w{Xm)jixi)jiXm) 



dxidxm JJ dPg{xk) 



+ 



- 1)2 



N-2 f w{xi)w{x2)jixi)jix2) 



N-l 



1 f w{xi)'^ 



dxi dx2 Y[ ^Pgi^k) 



k=3 



dPi^-i(x). 



N - 1 J^N-i S2,(x) 
The limiting behaviour of the first term is studied in the following Lemma: 
Lemma 2 With the assumptions and notation of Lemma 1, and defining 

w{xi)iv{x2)i{xi)j{x2) 



N-l 



we have 



■d2;idx2 W dP^(xfc) , 



A:=3 



(4.20) 



lim = \- w{y)j{y)dy 
N^oo \I Jn 

Proof: We follow the lines of the proof of Lemma 1. Defining again A]\f{e) := {x € ft^~^, |S'Ar(x) — /| < 
e} and A'j^{e) := \ AN{e), with < e < 1/2, we derive the estimates 



Aiv(e) 



'W{xi)w{x2) w{xi)w{x2) 



/2 



dPf-Hx) < ^ 



We 

T 
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and 



w{xi)w{x2) w{xi)w{x2) 



where the first term vanishes in the hmit N oo due to (4.18) and the second one by (4.19) with m = 3. 

■ 

By a slight modification of the above Lemma, one obtains the hmit of the second term of (4.20), namely 

lim / ^ dP^H^) = i / w{yf dPM . 



Therefore, the limit as iV ^ oo of (4.20) is 



w{y)3{y)dy \ , (4.21) 



V-i V (^-l)5iv(x) j " N^^ \lJn 

where we used Lemma 2. 

Remark 1 Lemmas 1 and 2 were formulated for the case k = 1, however, all the calculations can be 
easily generalized for any (fixed) value ofk. 

4.3 Derivation of kinetic and hydrodynamic description 

Let us denote p~^{t, x) := p{t, x, 1) and p^{t, x) := p{t, x, —1). We define u{t, x), the continuous analogue 
of the local average velocity (2.2), by 

where 

Soy,p-]{t) := j^w{\x-z\){p+{t,z)+p-{t,z))dz = 0^ . 

We extend the definition of u{t, x) to the whole domain Q by setting u{t, x) = for x G So\p~^ , p~]{t) , see 
Remark 2. 

The kinetic equation for p is obtained by setting = 1 in the BBGKY- hierarchy (4.10) with the 
molecular chaos assumption (4.13) and passing to the limit — t- oo. Using Lemma 1 and formula (4.21), 
we obtain 

q^{x;x,v) = p{t,x,v)u{t,x) , and r^{x;x,v) = p{t,x,v)u{t,x)'^ . (4.23) 

Consequently, (4.10) reduces to the following system of two equations 

dtp^ + d,p+ = -[^0 + bil - uf]p+ + [^0 + b{l + uf]p- , (4.24) 
dtp~-d^p- = -[7o + 6(l + n)V + [7o + fe(l-ii)>+. (4.25) 
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Equivalently, defining the mass density g and flux j by 

Q{t,x) ■.= p^{t,x) +p~{t,x) , j{t,x) ■.= p^{t,x) -p~{t,x) 

the system can be written in the hydrodynamic description as 



dtQ + dxj 
dtj + dxQ 

u{t, x) 



0, 



2[jo + b{l + u^)]j + 4bgu , 



X i 5o[£>](t) 



with 



^^W(\X - z\)Q{t,Z)^Z ' 

0, xG5o[i?](t) 



S y — z) dz = 



(4.26) 

(4.27) 
(4.28) 

(4.29) 



Remark 2 Due to the assumption [A2'], we have p'^{t,x) = p {t,x) = on So[p^,p ]{t). Therefore, 
both q^'^{x; x, v) in (4.14) and r^'i(a;; x,v) in (4.15) are equal to zero by definition. Consequently, 

q^{x;x,v)= lim q^'^{x; x,v) = , and r^{x;x,v)= hm r^'^{x; x,v) = , 
N^oo TV— >-oo 

and formulas (4.23) remain valid irrespective of the particular choice of the value ofu{t,x). This justifies 
our extension of the definition of u{t,x) by setting it equal to zero for x S So[p'^ , p~]{t) . 



5 Existence of solutions to the kinetic system 

The main goal of this Section is to prove the following Theorem: 

Theorem 1 Let 70 > 0, 6 > and w satisfy the assumptions [Al] and [A2']. Then, for every T > 
and every nonnegative initial datum {pq,Pq) G L°°(J7) x L'^{Q) there exists a nonnegative solution to 
the kinetic formulation (4.24)-(4.25) in L°°{[0,T] x 17). This also establishes solutions g = p'^ 
j = p~^ — p" of the hydrodynamic formulation (4.27) -(4.29) with the corresponding initial condition. 

Proof of Theorem 1: The proof is carried out in three steps and is only sketched here, omitting details 
where the techniques are standard. For notational convenience, we will work both with the kinetic and 
hydrodynamic representation of the system and treat {p'^,p~) and {g,j) as synonyms, related by (4.26). 

Step 1. First, we consider a linearized version of (4.24)-(4.25), where we solve for p~^ and p~ given a 
prescribed u € -L°°([0,T] x Q) with |u| < 1. This constitutes a strictly hyperbolic system with unique 
mild, nonnegative solution in C([0, T]; L°°(0)), constructed by a standard fixed point iteration (see, for 
instance, [10], Section 7.3). The solution is given by the Duhamel formula 

p+(t,x) = p+(0,x-t)+ / Q+[p+,p-]{s,x + {s-t))ds, (5.1) 

Jo 

p~{t,x) = p-{{),x + t)+ / Q~[p+,p~](s,x-(s-t))ds, (5.2) 



14 



with Q^\p^ ,p~] = T[lo + b{l — lb [70 + 6(1 + Moreover, for any fixed T > 0, we have apriori 

boundedness of p'^ and p~ in L°°{[0,T] x O), depending only on the initial condition. Indeed, denoting 
h(t) := sup^gQp+(i, x) + sup^^Q p~{t,x), and remembering the uniform boundedness \u\ < 1, we have 

h{t) < h{0) + (70 + 46) / h{s) ds , 

Jo 

and the apriori boundedness follows from an application of the Gronwall lemma on the time interval [0, T]. 
Step 2. We consider a regularized version of (4.24)-(4.25) where u is substituted by Ug, defined by 

J{x) 



Ue{x) :-- 



£ + R{x) 



with J{x):= j w{\x — z\)j{z) dz , R{x) := I w{\x — z\)q{z) dz . (5.3) 
Jo. Jo. 



For any fixed e > 0, a solution is found by the Schauder fixed point iteration on the mean velocity [10]. 
The compactness of the corresponding Schauder operator is provided by the Arzela-Ascoli theorem. In- 
deed, let us take a sequence vP with ||'w"||j;,oo([o T]xf7) — ^ ^'^^ P^'"^ ^"^^ P~'^ be the corresponding 
mild solutions of the kinetic system, given by the Duhamel formula (5.1)~(5.2) with in place of u, 
and let = p"*"'" + and j" = p^'" — As explained in Step 1, for any fixed T > we have 

IIp^'^II^cxiqo T]xn)i IIp^'"IIl°°{[o T]xn) bounded uniformly with respect to n. Defining the function 



U!{x) := / \w{z) — w{\z — x\)\dz for x G 0. , 
Jn 

one has lima;_^o ^^(2;) = (continuity of translation [21]). By Holder inequality, 

\J^ix)-J^iy)\<\\rh^^^^u;{x-y), 
with J'^{x) := j^w{\x — z|)j"(z)dz. Moreover, we have the uniform boundedness 



iL°°(n) 



\W\\t.i 



and analogous estimates hold for R^{x) := J^w{\x — z\)Q{z)dz. Consequently, we have 

< ^11/ 

e 



This uniform equicontinuity together with the uniform boundedness < 1 allows us to apply the 
Arzela-Ascoli Lemma and obtain the compactness of the Schauder operator. Therefore, for every fixed 
e > 0, we have a nonnegative solution {pf,pj) of the regularized system (4.24), (4.25), (5.3), uniformly 
bounded (with respect to e) in L°°([0,T] x il). 

Step 3. Finally, we pass to the limit e — > 0. Due to the uniform boundedness, a subsequence of pf p^ 
weakly* in L°°([0, T]; L°°(ri)); we need to show that the limit of the nonlinear terms QgU^ and u^je is gu 
and, resp., u^j, with u given by (4.22), or, equivalently, (4.29). The limit passage in the distributional 
formulation of the term with a test function (p G C^([0, T) x fi) is performed as follows: 



roo r rco r roo r 

/ / {QeUe - gu)Lpdxdt < / {Qe - g)uLp dx dt + 

JO Jn Jo Jn Jo Jn 



ge{ue — u)ipdxdt 
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The first term vanishes in the hmit e ^ due to the weak* convergence of Qs towards g with mp a vahd 
test function. Concerning the second term, we will show that Qeiue — u)ipdx tends to zero for almost 
all t € [0, T] and conclude the convergence of the time integral by the Lebesgue dominated convergence 
theorem. Let us fix t G [0, T] and define the set Ss by 

Ss := {x G n; R{x) < 5} , 

with (5 > and R defined in (5.3). Then, we have 



/ Qeiue - u)(pdx < 2 / 
JSs JS 



< 2 I ^e|v:'|dx > 2 / Q\ip\ 



dx 



= 2 / Q\ip\ dx , 

Jss\So <5->o 

where the second line is due to ^? = on 5o and because meas(55 \So) tends to zero as 5 — )■ 0. Next, for 
X e ^\Ss, 



R{x)Re{x) 



R{x) 



< -{\Re{x)-R{x)\ + \Je{x)-J{x)\), 

where Je(x) = J^w{\x — z\)ji;{z)dz and i?e(x) = f^w{\x — z\)qs{z) dz. Therefore, due to the uniform 
convergence of i?e and to R and, resp., J onil. (implied by the uniform equicontinuity and boundedness 
of the families {Re}e>o and {Je}e>o), we have 



/ 



QsiUs - u)ipdx 



< I supf^\5jn, - m 



Qe\^\ dx 



~ ^ v'-^^ ~ 11"'^ ~ ll'^'lli'Mf^) 



We conclude by passing 5^0. The limit passage in the term u^je is performed similarly (note that 

\je\ < Qe)- 



It is worth noting that the assumption [Al] of Theorem 1 can be relaxed. In fact, we posed the requirement 
[Al] of boundedness of w on VL in order to establish the definition (2.2) of uf" in the discrete model and 
pass to the limit N ^ oo. However, at the level of the kinetic or hydrodynamic description, we may 
relax this to w (z L^(r2). Moreover, formally it is possible to consider even singular weights, in particular 
w = (5o, which leads to u = j/ q and removes the nonlocality. In fact, one can see the choice w = 5q as the 
limiting case when the interaction radius shrinks to zero: for almost all x G such that ^(x) 7^ 0, one has 



lim 



Lx[o,<7](|a;-^;|)j(2)d2; _ j(x) 



"^^o !n X[o,<7] {\x - z\)q{z) dz q{x) ' 

where X[o,(t] is the characteristic function of the interval [0,cj]. One can interpret this as a model where 
only pointwise local observations of the system are possible. 
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By a slight modification of the proof of Theorem 1 it is possible to show that given a sequence of 
weights Wn converging strongly in L^(il) to w, the solutions {pn,jn) corresponding to the weights Wn 
converge weakly* in L°°{[0,T]; L°°{il.)) to the solution corresponding to the weight w. However, we need 
the condition [A2'] to be satisfied uniformly; consequently, the question whether and how the solution 
corresponding to w = Sq can be obtained as a limit of solutions with Wa = X[o,a] as it — )• remains an 
interesting open problem. An even more interesting question is what is the limit of the discrete model as 
— 7> oo if the interaction radius is shrinking as some power of 1 /N. This question is studied in Section 7 
below. 

6 Long time behaviour 

In this Section we provide several conjectures about the long time behaviour of the kinetic system (4.22), 
(4.24), (4.25) or, equivalently, the hydrodynamic system (4.27)-(4.29). To get some insight into the long 
time dynamics, we start with a numerical example. We solve the kinetic system using standard semi- 
implicit finite difference method with upwinding. The initial condition is Pq = 2.2 on [0.125,0.375] and 
zero otherwise, Pq = 1.8 on [0.625,0.875] and zero otherwise. In Figure 6 we show the results for the 
choice of parameters 6 = 1 and 70 = 0.3 and the weight function w = X[o,o.2]- 

We conjecture that with regular weights w satisfying [Al] and [A2'], the solutions (p, j) to (4.27)-(4.29) 
converge to the constants p = 1 and j = jg with some jg G M, exponentially fast as t 00. Moreover, 
in the large noise case 70 > b, we hypothesize that jg = 0. Unfortunately, we are only able to provide an 
analytical proof in the rather special case w = 1 and 70 > b: 

Lemma 3 Assuming w = 1, we have u{t,x) = u{t), where u{t) satisfies the ordinary differential equation 

tt = -2(7o + 6(n2 - l))u, (6.1) 
subject to the initial condition u{0) = J^j(0,x)dx. Moreover, 

(i) if 70 < b, then Vrnit-, 00 u{t) = sign(n(0)) y^l - 7o6"i, 

(ii) if Jo > b, then limt^aDu{t) = and 

|n(t)| < |u(0)|e-2(To-'')*. (6.2) 
Moreover, j converges to zero exponentially fast in the L"^ -sense: 

[ f{t,x)dx < ce-^To* 

for a suitable constant c. 

Proof: Integrating (4.28) with w = 1 with respect to x € we obtain (6.1). This has the stationary 
state n = 0, which is stable if and only if 70 > b. Moreover, if 70 < b, two additional stable stationary 
states u = zby^l — 7o6^^ exist. This establishes the first statement. Further, we have 

^ Qn^) = -2(70 + biu' - 1)W < -2(70 - bW , 
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Figure 3: Numerical results with w = X[o,o.2]) b = 1 and 70 = 0.3: p+ and p converge to constant states, 
p converges to 1. 
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and an application of the Gronwall lemma gives (6.2). 

To prove the convergence of j to zero, we consider the identity 



1^ 
2 dt 



/ (e*^ + j^) dx = -2(70 + 6(1 + n^)) / f dx + ibu / gj dx . (6.3) 
Jn Jn Jn 

An application of the Cauchy-Schwartz inequality and the decay rate (6.2) yield 

Abu I gjdx<2b [ fdx + 2bu^ [ g"^ dx < 2b [ dx + 26|w(0)pe-^(^0"^)* [ g"^ dx . 
Jn Jn Jn Jn Jn 

Then, an integration of (6.3) in time leads to 

[ ^2(r,x)dx < co + 46|^/(0)pe-^(^«~^)* /" [ g'^it,x)dxdt , 
Jn Jo Jn 

with Co := gQ{x) + jQ{x) dx. Consequently, by the Gronwall lemma, g^ dx is bounded uniformly in 
time by a constant ci if 70 > b. Inserting this information into (6.3) gives 

[ f{T,x)dx < co-470/" [ f{t,x)dxdt + ibci\u{0)\'^ [ 6-^(^0-'')* dt 
Jn Jo Jn Jo 

< c-470 /" [ f{t,x)dxdt 
Jo Jn 

and an application of the Gronwall lemma yields the second statement. 



6.1 The case w = 60 

In this subsection we briefly discuss the long time behaviour in the singular case w = 5o. Again, we start 
with a numerical example where we solve the kinetic system with the parameters b = 1 and 70 = 0.3. The 
initial condition is chosen as before, see Figure 6 (left panels). In Figure 6.1 we present the time evolution 
of p~^, p~ p, j and u. Based on the numerical observations, we conjecture that, for the small noise case 
7o < b, the long time dynamics are given by the travelling wave profiles {pf,pj), satisfying 

dtpt + dxPt = , 
dtP'^ - dxPj = , 



and ^1 G {uc,, —Ug} with Us = \/l — 706 Then, it is a matter of a simple consideration to deduce 

Ps ~^Ps 

that one of the functions, say pj , has to be a global constant, while the other one, pf, is a piecewise 
constant assuming only two values {p^i,pt2}^ satisfying the relations 



P~^ / 1 w \ ^ 
Pt,pt2 = iPsf and = 



These relations chracterize the dynamic equilibria between the densities of individuals marching to the 
left and to the right, in dependence on the parameter values 70 and b. 
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Figure 4: Numerical results with the singular weight w = 6o, b = 1 and 70 = 0.3: p converges to a 
constant, while p'^ becomes a piecewise constant travelling wave profile and u jumps between the values 
±Us with Us = \l\ - 706-1 ^ 0.8367. 
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For the large noise case (70 > b), we prove that the asymptotic state is again i = and p = 1. In 
fact, j converges to in the L^-sense exponentially fast as t — > 00. This follows from an application of 
the Gronwall lemma to the estimate 

< -2(70-6) j fdx. 

7 The effect of shrinking interaction radius 

In biological applications, the weight function w is often chosen as the characteristic function of the 
interval [0, a], where a > is the interaction radius. It is important to study the dependence of collective 
behaviour on the size of o". In this section, we will consider a theoretical limit o" — )■ 0. Clearly, letting a ^ 
with fixed N leads to a trivial model without any interactions (it corresponds to the choice b = 0). On 
the other hand, letting first N — ?• 00, followed by o" ^ 0, we obtain the model with w = 6q, as mentioned 
in Section 5. Consequently, the limits o" ^ and — )• 00 do not commute. From the modelling point 
of view, it is natural to study the limit N ^ 00 with the interaction radius shrinking as A^""^ for some 
a > 0. Indeed, as the space is getting more crowded, the visibility is reduced and every individual can 
only take into account its closest neighbours. Actually, taking w = X[o,N-'^] letting — )■ 00, we will 
show that a significant limit is obtained with a = 1, leading to a new kinetic model. The results are 
summarized in the following Theorem. 

Theorem 2 Let w = X[o,n~"]- Then, the formal limit as N ^ 00 of the BBGKY hierarchy (4.10) is 
(1) For a > 1: 

d^Q + 2{jo + b)dtQ = d^Q. (7.1) 



(2) ForO<a< 1: 



(3) For a = l: 



dtQ + d.^3 = 0, (7.2) 
dtj + d^Q = -2(^7o + 6(^l + ^))i + 46j. (7.3) 



dtQ + d^j = 0, (7.4) 
dtQ-d.,j = -2(7o + 6(l + r?))j + 46j(l-exp(-^)), (7.5) 



with 



7] = — {I- exp(-^)) + (1-^1 exjp{-g) [Et{Q) - 7 - ln{Q)] , (7.6) 

where Ei(z) = fl^ exp(— s) ds is the so-called exponential integral function and 7 = 0.577 ... is 
Euler's constant. 
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We see that for < a < 1, we obtained the kinetic model with u = j/p (i.e., the same as if one would 
first pass to — t- oo and then o" — )> 0). If q > 1, we get the model with no interactions (i.e., as if 
one would first pass to (j — )> and then — )> oo) which is described by the telegraph equation for g 
[19]. In the significant limit with a = 1 we obtained a new model, which is in fact the hydrodynamic 
model (4.27)-(4.28) with the function ^(1 — exp(^)) in place of u and with rj, given by (7.6), in place of 



Proof: All we need to do is to recalculate the limits A'" — )• cx) in the expressions (4.11) and (4.12) with 
w = X[o,N-<^]- Assume that x = is a Lebesgue point of g{x) and j{x), i.e., 

lim- r g{y)dy = g{0), (7.7) 



cr-S-O a JQ 

and similarly for j. Moreover, assume that g{0) > 0. In the same way as in Lemma 1, we calculate 

^ \ ^ / 2^m=lXlO,N-°']{Xm)Vm yj , x , 

Qn := 2^ / — ^lv3i -— [[p{xi,Vi)d^ 

V6{±1}^-1-^^ Em=lX[0,7V-](2;m) 
2^rn=lXlO,N'"]{Xm) 

= {N-l)Kr,[ , ^ -- dP^-\^), 

JnN-2 1 + XlO,N-"]i.Xm) 

where we denoted 



IN 



Obviously, the sum X]m=i X[o,N-'^]ixm) only takes the values A; = 0, . . . , A^ — 1, and 

N) 



Therefore, 



k=0 

1 



1 + ^^^^ XlO,N—]iXm) 

Due to (7.7), we have limAr_).oo N°'Ij\f = g{0) and limjv_^oo N'^Kjy = Therefore, 

hm Qm= hm (1-(1- 7^)^-1) =M (l_(l_/^)A^-i) , 

and 



lim (1 - (1 - In) ) = { 1 for < a < 1 , 

N—>-oo ' 



1 - e"^'(°) for a = 1 , 
1 for < a < 

for a > 1 . 
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Using the same technique as above, one calculates that 



„ X-- f fj2m=lX[0,N-^]iXm)Vm\ y-r 

vg{±l}iv-i -^f^^-' V Em=lX[0,N-'^]iXm) J 

V -'TV / 



k=l ^ ^ 

It is easily shown that Gn{In) vanishes in the limit N ^ oo whenever a ^ 1, while with a = 1 the law 
of rare events ([11]) gives 



with 



hm Gn{In) = e-^(°) E %f = [EKt'CO)) - 7 " HQiO))] , 

k=l 

We assume that q and j be integrable functions, so that almost every x G O is a Lebesgue point for them. 



Remark 3 Our analysis can he seen as a first step towards a so-called topological model, where the agent 
interactions are based on some connectivity graphs. For instance, one can consider the situation where 
every agent interacts only with its nearest neighbour. In this case, we have a different time dependent 
weight function Wi for every agent, namely Wi = X[o,ai] with ai = miumj^i \xi — Xm\- Then, the passage to 
the corresponding continuum model is a completely open problem; however, let us observe that 



— 1 ^ \xi — a 



-i/p 



mm \ Xi 



N — I ^ \xi — Xm\^ j p-^oo m 

Therefore, it would be interesting to consider the discrete system with Wi = X[o,o-i]' 

and study the limit N ^ oo and p ^ oo (possibly with p = N). Moreover, let us observe that "on 
average", \xi — Xm\ ~ j^- Therefore, 

-i/p 



so we are in the situation of the significant limit a = 1 described above, and we might believe that the new 
kinetic model obtained in this limit can have some connection with the limit N ^ oo and p — ?■ oo o/ (7.8). 
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8 Discussion 



We introduced an individual based model with velocity jumps aimed at explaining the experimentally 
observed collective motion of locusts marching in a ring shaped arena [2]. The frequency of individual 
velocity jumps increases with a local or global loss of group alignment. We showed that our model has the 
same predictive power as the model of Czirok and Vicsek, in particular, it exhibits the rapid transition to 
highly aligned collective motion as the size of the group grows and the switching of the group direction, 
with frequency rapidly decreasing with increasing group size. Moreover, in the limit — > oo we obtained 
a system of two kinetic equations. We proved existence of its solutions and a partial result about the 
long time behaviour. Finally, we studied the effect of shrinking the interaction radius a in the discrete 
model as the number of individuals, A^, tends to infinity. We showed that in the significant limit where 
a shrinks as one obtains a new kinetic model. 

Kinetic approach has previously been used in the literature to understand collective dynamics of 
individual based models. Carrillo et al [3] found the double milling phenomena in the kinetic formulation 
of the model of self propeled particles with three zones of interactions. The kinetic description of the 
Cucker-Smale model was introduced in [14], which can also be be derived from the Boltzmann-type 
equation, see [4], or Povzner-type equation, [12]. For the survey of the most recent results see the 
review [5]. 

Several interesting questions remain open, offering space for future investigations. For example, the 
kinetic system (4.24)-(4.25) deserves a better analysis, in particular, uniqueness of solutions and more 
complete investigation of the long time behaviour. It would also be interesting to know if and how the 
solutions corresponding to w = Sq can be derived as a limit of solutions corresponding to w = X[o,a] 
(T — 7> 0. Another interesting direction of future research was formulated in Remark 3. 
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